!
! statistics -  a library of basic statistical functions
!

module statistics

  implicit none

  integer, parameter, private :: dp = selected_real_kind(15)

  private :: rm, qmed

contains


  ! arithmetical mean, standard deviation
  subroutine mean(x,t,dt)

    implicit none
    real(dp), intent(in) :: x(:)
    real(dp), intent(out) :: t,dt
    integer :: n

    n = size(x)
    if( n == 0 ) then
       t = 0.0_dp
       dt = 0.0_dp
    else if( n == 1 ) then
       t = x(1)
       dt = 0.0_dp
    else
       t = sum(x)/n
       dt = sqrt(sum((x-t)**2)/(n - 1))
    end if

  end subroutine mean


  ! median
  function median(x) result(t)

    real(dp), dimension(:), intent(in) :: x
    real(dp), dimension(size(x)) :: y
    real(dp) :: t
    integer :: n

    n = size(x)
    y = x
    call qmed(n,y,n/2+1,t)

  end function median

  ! robust mean, standard deviation
  subroutine rmean(x,t,dt)

    use cutoff

    real(dp), dimension(:), intent(in) :: x
    real(dp), intent(out) :: t,dt

    call rm(x,t,dt,tukey,dtukey,50,0.1_dp,sqrt(epsilon(t)))

  end subroutine rmean

  subroutine rm(x,t,dt,psi,dpsi,maxit,relerr,abserr)

    interface
       function psi(x)
         integer, parameter :: dp = selected_real_kind(15)
         real(dp) :: psi,x
       end function psi
       function dpsi(x)
         integer, parameter :: dp = selected_real_kind(15)
         real(dp) :: dpsi,x
       end function dpsi
    end interface

    integer, optional, intent(in) :: maxit
    real(dp), optional, intent(in) :: relerr, abserr
    real(dp), dimension(:), intent(in) :: x
    real(dp), intent(out) :: t,dt

    real(dp), dimension(size(x)) :: xx
    integer :: n,i,it
    real(dp) :: d,r,s,s2,sum1,sum2,sum3,rp

    n = size(x)
    if( n < 2 ) return

    ! initial estimation of parameters
    xx = x
    t = median(xx)
    xx = abs(x - t)
    s = median(xx)/0.6745
    s2 = s**2*n/(n-1)

    ! identical points on input
    if( abs(s) < epsilon(s) )then
       call mean(x,t,dt)
       return
    endif

    it = 0
    do
       it = it + 1

       ! corrector's estimation 
       sum1 = 0.0 
       sum2 = 0.0
       sum3 = 0.0
       do i = 1,n
          r = (x(i) - t)/s
          rp = psi(r)
          sum1 = sum1 + rp
          sum2 = sum2 + dpsi(r)
          sum3 = sum3 + rp**2
       enddo

       if( abs(sum2) < epsilon(sum2) ) exit

       d = s*sum1/sum2
       t = t + d
       dt = s2*sum3/sum2**2

       ! exit of iterations:
       if( present(relerr) .and. d**2/dt < relerr**2 ) exit
       if( present(abserr) .and. abs(d) < abserr ) exit
       if( present(maxit)  .and. it > maxit ) exit
       if( .not. present(maxit) .and. it > 100 ) exit
       
    enddo

    ! deviation's estimation
    dt = sqrt(dt*n)

  end subroutine rm

  subroutine qmed(n,a,k,x)

    integer, intent(in) :: n,k
    real(dp),intent(in out) :: a(:)
    real(dp),intent(out) :: x
    real(dp) :: w
    integer :: l,r,i,j

    l = 1
    r = n
    do while( l < r )
       x = a(k)
       i = l
       j = r
       do
          do while( a(i) < x )
             i = i + 1
          enddo
          do while( x < a(j) )
             j = j - 1
          enddo
          if( i <= j ) then
             w = a(i)
             a(i) = a(j)
             a(j) = w
             i = i + 1
             j = j - 1
          endif
          if( i > j ) exit
       enddo
       if( j < k ) l = i
       if( k < i ) r = j
    enddo

  end subroutine qmed

end module statistics
